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(Received  18  March  2010;  published  12  luly  2010) 

Biological  membranes  are  involved  in  numerous  intriguing  biophysical  and  biological  cellular  phenomena  of 
different  length  scales,  ranging  from  nanoscale  raft  formation,  vesiculation,  to  microscale  shape  transforma¬ 
tions.  With  extended  length  and  time  scales  as  compared  to  atomistic  simulations,  solvent-free  coarse-grained 
membrane  models  have  been  exploited  in  mesoscopic  membrane  simulations.  In  this  study,  we  present  a 
one-particle-thick  fluid  membrane  model,  where  each  particle  represents  a  cluster  of  lipid  molecules.  The 
model  features  an  anisotropic  interparticle  pair  potential  with  the  interaction  strength  weighed  by  the  relative 
particle  orientations.  With  the  anisotropic  pair  potential,  particles  can  robustly  self-assemble  into  fluid  mem¬ 
branes  with  experimentally  relevant  bending  rigidity.  Despite  its  simple  mathematical  form,  the  model  is  highly 
tunable.  Three  potential  parameters  separately  and  effectively  control  diffusivity,  bending  rigidity,  and  spon¬ 
taneous  curvature  of  the  model  membrane.  As  demonstrated  by  selected  examples,  our  model  can  naturally 
simulate  dynamics  of  phase  separation  in  multicomponent  membranes  and  the  topological  change  of  fluid 
vesicles. 

DOI:  10.11 03/Phy sRe vE. 82.0 11905  PACS  number(s):  82.70.Uv,  87.17.Aa,  81.16.Dn 


I.  INTRODUCTION 

Due  to  the  hydrophobic  interactions,  amphiphilic  mol¬ 
ecules  in  aqueous  solutions  self-assemble  into  monolayer 
(e.g.,  amphiphilic  block  copolymers  [1])  or  bilayer  mem¬ 
branes,  where  rodlike  molecules  align  their  longitudinal  axes 
parallel  to  each  other  and  are  able  to  diffuse  laterally  at 
physiologically  relevant  temperatures.  While  atomistic  simu¬ 
lations  [2]  have  been  routinely  performed  to  elucidate  physi¬ 
cal  mechanisms  involved  in  localized  lipid  organizations, 
they  are  inaccessible  to  many  intriguing  membrane  processes 
[3,4]  that  occur  on  the  length  scale  much  larger  than  mem¬ 
brane  thickness.  To  overcome  length  and  time  scale  limita¬ 
tions  of  atomistic  simulations,  significant  efforts  have  been 
devoted  to  develop  particle-based  coarse-grained  models, 
where  each  lipid  molecule  is  coarse-grained  into  several  con¬ 
nected  beads.  In  general,  particle-based  models  can  simulate 
membrane  fluidity,  hydrodynamic  effect,  and  membrane  to¬ 
pological  changes  naturally. 

In  the  explicit-solvent  coarse-grained  models  [5-7],  sol¬ 
vent  molecules  are  also  coarse  grained,  i.e.,  grouping  several 
water  molecules  into  a  coarse  grain.  Explicit-solvent  coarse¬ 
grained  models  stabilize  membrane  in  fluid  phase  by  explic¬ 
itly  accounting  for  the  hydrophobic  interactions  between  sol¬ 
vent  and  lipid  particles.  While  explicit-solvent  strategy  is 
convenient  and  natural,  yet  it  comes  with  a  high  computa¬ 
tional  cost.  In  coarse-grained  simulation  settings,  membrane 
particles  spread  on  a  two-dimensional  (2D)  surface,  whereas 
solvent  particles  occupy  the  bulk.  The  solvent  degrees  of 
freedom  thus  vastly  outnumber  the  lipid  particles  even  for  a 
modest  sized  membrane,  limiting  the  accessible  length  and 
time  scales  of  explicit-solvent  coarse-grained  models. 
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With  the  limitations  of  the  explicit-solvent  models  in 
mind  solvent-free  models  are  highly  desired.  For  solvent-free 
models  in  which  water  molecules  are  not  explicitly  modeled, 
effective  intermolecular  attractive  interactions  need  to  be  de¬ 
vised  to  substitute  for  the  hydrophobic  interactions  between 
water  molecules  and  hydrocarbon  chains  of  amphiphilic  mol¬ 
ecules,  which  constitutes  the  major  challenge  in  the  develop¬ 
ment  of  solvent-free  membrane  models.  Such  intermolecular 
interactions  have  been  successfully  developed  in  chain-of- 
bead  models,  where  each  amphiphilic  molecule  is  coarse 
grained  into  a  chain  of  multiple  connected  beads  [8-13]. 
These  chain-of-bead  models,  e.g.,  the  three -bead  model  by 
Cooke  and  Deserno  [11,14,15],  have  been  exploited  in  me¬ 
soscopic  membrane  simulations. 

To  push  the  length  and  time  scales  to  the  highest  possible 
level  while  retaining  key  membrane  properties,  particle- 
based  membrane  models,  where  a  single  or  a  cluster  of  am¬ 
phiphilic  molecules  are  represented  by  only  one  particle,  are 
highly  desired.  However,  attempts  to  develop  such  models 
are  often  frustrated  by  the  complexity  of  mimicking  the  hy¬ 
drophobic  effects  in  the  absence  of  explicit  solvents  [4].  In 
the  pioneering  model  of  Drouffe  et  al.  [16],  particles  are  able 
to  self-assemble  into  one -particle-thick  fluid  phase.  How¬ 
ever,  membrane  bending  rigidity  of  this  model  is  signifi¬ 
cantly  lower  than  the  experimental  range.  In  addition,  its  use 
of  multibody  potential  complicates  the  model  and  increases 
the  computational  cost.  Along  this  direction,  several  other 
models  were  proposed,  aiming  to  avoid  using  multibody  po¬ 
tentials  or  to  reproduce  experimentally  relevant  membrane 
properties  [17-19].  In  a  recent  model  of  such  type,  Kohyama 
[19]  extended  the  model  of  Drouffe  el  al.  [16]  to  a  pair 
potential,  where  the  bending  rigidity  is  controlled  by  a  time- 
dependent  variable.  Despite  varying  degrees  of  success  of 
these  models,  it  remains  unclear  whether  an  anisotropic  pair 
potential  as  simple  as  Lennard-Jones  (LJ)  potential  for  the 
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FIG.  1.  (Color  online)  Schematics  of  the  orientation-dependent 
interparticle  interaction.  Kinetically,  each  particle  is  axisymmetric 
with  a  particle-fixed  unit  vector  n  representing  the  axis  of  symme¬ 
try,  and  with  mass  of  m.  The  interparticle  interaction  is  both  dis¬ 
tance  and  orientation  dependent.  The  angle  90  is  a  model  parameter 
charactering  the  spontaneous  curvature.  The  configuration  corre¬ 
sponding  to  9j=6j=8a  is  made  to  be  the  most  energetically  favor¬ 
able  relative  orientation  between  two  particles.  The  two  halves  of 
each  particle  are  colored  distinctly  to  indicate  the  orientation  of  the 
particle.  The  spontaneous  curvature  c0  is  related  to  60  via  c0 
~2  sin  dgldo,  where  d0  is  the  average  interparticle  distance. 

isotropic  case  exists  to  mimic  the  hydrophobic  effect  while 
yielding  biologically  relevant  properties.  In  addition,  for 
large  length  scales,  triangulated  membrane  models  [20-23] 
have  been  intensively  used  in  simulating  biological  mem¬ 
branes  and  protein  networks  [24-26].  The  finite  element 
method  [27],  the  meshless  method  [28],  and  the  phase  field 
approach  [29-31]  have  also  been  applied  to  study  fluid  mem¬ 
branes. 

The  present  work  highlights  a  one-particle-thick  fluid 
membrane  model.  The  interparticle  interaction  is  described 
by  a  soft-core  pairwise  potential,  where  the  interaction 
strength  is  dependent  on  the  relative  orientations.  The  pair¬ 
wise  interparticle  interaction  potential  is  very  simple  in  its 
mathematical  form  but  highly  tunable,  able  to  robustly  self- 
assemble  into  fluid  membranes  with  biophysically  relevant 
properties.  From  length  and  time  scale  mappings,  each  par¬ 
ticle  in  our  model  represents  a  few  lipid  molecules  in  the 
lateral  direction.  Due  to  its  high  computational  efficiency, 
our  model  enables  simulations  of  large-scale  membrane  to¬ 
pological  changes  and  phase-separation  dynamics,  as  demon¬ 
strated  by  the  selected  examples. 

II.  COARSE-GRAINED  MODEL 

The  coarse-grained  particles  in  our  model  are  axisymmet¬ 
ric,  with  their  axes  of  symmetry  representing  the  longitudinal 
direction  of  lipid  molecules.  The  model  features  a  soft-core 
pairwise  interparticle  potential  with  the  interaction  strength 
weighed  by  the  relative  orientations  of  the  particle  pair.  Cor¬ 
respondingly,  the  interparticle  potential  is  constituted  of  two 
functions,  u(r )  and  (/>(f,;, n,,^),  which,  respectively,  de¬ 
scribe  the  distance  and  orientation  dependences.  Figure  1 
depicts  a  generic  relative  position  and  orientation  of  such  a 
particle  pair,  where  the  two  halves  of  each  particle  are  col¬ 
ored  distinctly  to  indicate  its  orientation.  This  coloring 


scheme  facilitates  visualization  [32]  in  simulations  presented 
later.  We  denote  r,  and  r  .■  the  center  position  vectors  of  par¬ 
ticles  i  and  j,  respectively.  The  interparticle  distance  vector  is 
then  rij=ri-rj.  We  also  denote  r=  |r,;|  and  r;/=ri;/ r.  The  unit 
vectors  n,  and  n;  represent  the  axes  of  symmetry  of  particles 
i  and  j,  respectively.  For  simplicity,  the  rotational  degree  of 
freedom  about  the  axis  of  symmetry  of  each  particle  is  ne¬ 
glected.  Therefore,  each  particle  carries  five  degrees  of  free¬ 
dom,  i.e.,  three  translational  and  two  rotational  degrees  of 
freedom  (noting  the  constraint  n-n=l). 

In  searching  for  the  functional  form  of  the  distance- 
dependent  function  u(r),  we  found  that  the  classical  12-6  LJ 
potential  only  leads  to  solid  membranes  at  low  temperatures 
and  gas  phase  at  high  temperatures.  Missing  of  a  fluid  phase 
in  between  is  due  to  the  steep  energy  landscape  of  the  12-6 
LJ  potential:  the  interparticle  interaction  forces  are  either  too 
strong  (near  the  equilibrium  distance)  to  permit  particle  dif¬ 
fusion  or  too  weak  (once  escaped  from  the  equilibrium  dis¬ 
tance)  to  hold  particles  together.  In  order  to  be  able  to  tune 
the  restoring  force,  we  adopt  the  following  two-branch  func¬ 
tion  «(r): 


u(r)  = 


-  e  cos2f 


TT  {r~  r miJ 

.2  (rc  -  r min)  J  ’ 


r  <  r 


min 


<r<rc 


(1) 


where  e  and  a  are  the  energy  and  length  units,  respectively. 
The  repulsive  branch  uR{r)  adopts  the  4-2  LJ  type  potential 
that  has  a  much  lower  restoring  force  than  the  12-6  LJ  po¬ 
tential  at  the  equilibrium  distance.  The  attractive  branch 
uA{r)  is  of  a  cosine  function  that  smoothly  decays  to  zero  at 
the  cutoff  radius  rc.  The  exponent  £  tunes  the  slope  of  the 
attractive  branch  [see  Fig.  2(a)]  and  hence  the  diffusivity  of 
the  particles.  The  two  branches  smoothly  meet  at  r=rmin  with 
C1  continuity.  We  set  the  distance  at  the  minimum  of  the 
potential  rlnln=  v2cr,  the  same  as  that  in  the  12-6  LJ  poten¬ 
tial,  and  rc=2.6cr  to  include  second-neighbor  interactions. 

Following  the  treatment  in  anisotropic  potentials  for  liq¬ 
uid  crystal  or  colloids  [33,34],  we  use  4>  to  weigh  the  inter¬ 
action  strength  for  different  relative  orientations,  leading  to 
the  final  form  of  the  anisotropic  pair  potential: 

rr,  .  I  uR(r)  +  [l  ~  ^(ri-,n,-,n-)]e,  r  <  rmm 

f/(rl7,n;,n)  =  i 

[wA(r)0(fi7,n,,n,),  rmin  <  r  <  rc 

(2) 


The  orientation-dependent  function  scales  the  attractive 
branch  of  the  distance-dependent  potential,  while  shifts  up¬ 
ward  the  repulsive  branch  [shown  in  Fig.  2(b)].  The  separate 
operations  ensure  a  fixed  distance  of  minimum  energy  for 
different  relative  orientations  between  two  particles.  The 
orientation-dependent  function  substitutes  for  the  hydropho¬ 
bic  effects  and  takes  the  following  form: 


4>=  1  +  /r[fl(r,7,n,,n7)  -  1],  (3) 

where 
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FIG.  2.  (Color  online)  (a)  The  slope  of  the  attractive  branch  of  the  distance-dependent  function  varies  with  the  exponent  £.  The  solid 
black  curve  represents  the  repulsive  branch  of  the  function,  (b)  The  interparticle  potential  U  as  a  function  of  r  for  three  different  6t  (6j 
=  0j)  with  £=4,  /x  =  3,  and  dQ=0.  The  double  arrows  denote  the  orientations  of  a  particle  pair. 


a  =  (n,-  X  fy)  •  (ny-  X  ri;)  +  sin  0o(n j  -  n,)  •  ri;-  -  sin2  0O. 

(4) 


constraint  n,-n,=  l,  the  governing  equations  for  n,-  can  be 
derived  using  the  Lagrange  equations  with  constraint  forces. 


The  function  a  (hence  (f> )  reaches  its  maximum  of  1  when 
dj=6j=dQ  (see  Fig.  1  for  0h  Q-r  and  0O),  and  is  less  than  1 
otherwise.  Hence,  the  relative  orientation  0,  =  ();  =  0O  is  most 
energetically  favored.  Since  90  specifies  the  favorable  inter¬ 
particle  orientations,  it  directly  links  to  the  spontaneous  cur¬ 
vature  of  the  model  membrane.  The  parameter  /x  weighs  the 
energy  penalty  when  the  particles  are  disoriented  from  90 , 
and  is  thus  related  to  the  bending  rigidity  of  the  model  mem¬ 
brane. 

Corresponding  to  the  two  sets  of  degrees  of  freedom,  par¬ 
ticle  center  positions  and  orientations,  there  are  two  sets  of 
the  equations  of  motion  for  the  coarse-grained  model.  The 
first  set  governs  the  time  evolution  of  the  particle  center  po¬ 
sitions. 


where  m ,■  is  the  mass  of  particle  i,  Ui=1,jU{rij ,n;,n7),  and  j 
runs  over  all  the  neighbors  of  i.  The  second  set  of  equations 
governs  the  time  evolution  of  the  particle  orientation,  which 
can  be  derived  from  Euler’s  rigid  body  dynamics  equations. 
However,  considering  that  there  are  only  five  degrees  of  free¬ 
dom,  the  equations  of  motion  governing  particle  orientations 
can  be  derived  in  a  more  efficient  manner.  We  treat  n,  of 
particle  i  as  three  generalized  coordinates  with  a  geometric 


OUj 

Ifli  =  -  —  +  Mi,  (6) 

r)ni 

where  /,  is  the  moment  of  inertia  (/,  is  fixed  to  1  •  nijCr  in  this 
work),  A.,-  is  the  Lagrange  multiplier  and  has  the  following 
relation  with  n,-  and  ri„ 

dUj 

h=  —  •  n/ -/A'  ii,,  (7) 

fin, 

Our  coarse-grained  molecular  dynamics  (CGMD)  simula¬ 
tions  presented  below  for  planar  membranes  are  performed 
in  the  N2,T  ensemble,  where  X  is  the  membrane  tension. 
Simulations  for  vesicles  are  performed  in  the  NVT  ensemble. 
We  adopt  the  Nose-Hoover  thermostat  [35,36]  to  maintain 
the  system  at  desired  temperatures.  Since  our  model  is 
solvent-free,  the  rigid-body  translational  and  rotational  mo¬ 
tions  are  removed  at  each  time  step  in  our  simulations,  which 
may  otherwise  cause  significant  numerical  errors.  For  a  pla¬ 
nar  membrane  with  periodic  boundary  conditions,  the  Ber- 
endsen  pressure  coupling  algorithm  [37]  is  modified  to  a  2D 
case  to  maintain  a  desired  tension  Xo  in  the  membrane  by 
rescaling  the  particle  coordinates  and  box  size  at  each  time 
step.  The  scaling  factor  \  is 


FIG.  3.  (Color  online)  Snapshots  from  CGMD  simulations  demonstrating  self-assembly  of  randomly  distributed  particles  into  vesicles. 
The  parameters  used  in  the  simulations  are  £=4,  /U  =  3,  0a  =  0,  and  kBr=0.17e. 
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FIG.  4.  (Color  online)  Phase  diagram  in  the  (£,T)  plane  at  zero 
tension.  Three  regions  representing  gel,  fluid,  and  gas  phases  are 
identified,  separated  by  solid  lines.  A  broad  fluid  phase  region  with 
diffusion  constant  on  the  order  of  O.Ict2/ t  exists. 

x=  1  +  r-^r[So-2(0],  (8) 

3  tpKa 

where  rp  is  the  relaxation  time,  KA  is  membrane  area  com¬ 
pression  modulus,  At  is  the  time  step. 

III.  MEMBRANE  PROPERTIES 

We  next  demonstrate  that  our  interparticle  potential  cap¬ 
tures  well  the  hydrophobic  interactions  and  can  lead  to  ro¬ 
bust  self-assembly  of  particles  into  fluid  membranes.  Starting 
with  a  random  distribution  of  5882  particles  enclosed  in  a 
cubic  box  with  periodic  boundary  conditions,  CGMD  simu¬ 
lations  were  performed.  The  particle  system  is  thermostated 
at  a  constant  temperature  of  kHk=().  I 7k,  with  kH  the  Boltz¬ 
mann  constant  and  T  the  temperature.  Figure  3  depicts  four 
simulation  snapshots.  At  intermediate  stages  (r=2000r, 
where  t=  a\ml s  is  the  time  scale),  membrane  flakes  are 
formed.  The  membrane  flakes  then  coalesce  into  large  planar 
membranes  (r=5000r),  and  finally  close  to  form  vesicles  due 
to  the  edge  effect  (f=12000r).  In  our  simulations,  the  soft¬ 
core  interaction  potential  allows  us  to  adopt  a  relatively  large 
time  step  (Ar=0.02r). 


We  first  characterize  the  diffusion  constant  as  a  function 
of  the  exponent  £  and  temperature  T.  The  preassembled 
square  membrane  adopted  for  this  study  consists  of  5822 
particles  with  side  length  ~70cr.  We  maintain  zero  mem¬ 
brane  tension  in  our  simulations  and  set  /i  =  3  and  0Q=0.  The 
in-plane  diffusion  constant  D=(s])l4t  is  systematically  cal¬ 
culated  at  zero  membrane  tension,  where  s ,■  is  the  diffusion 
distance  of  particle  i  over  time  period  t.  This  allows  us  to 
construct  a  phase  diagram  of  the  diffusion  constant  on  the 
(£,  T)  plane,  as  shown  in  Fig.  4.  The  particle  membrane  is 
considered  in  gel  phase  if  D  <  0.01  a2!  r,  and  in  gas  phase  if 
at  least  one  particle  flies  away  from  the  membrane  during  the 
simulations.  We  identified  a  broad  fluid  phase  region  in 
which  the  diffusion  constant  is  on  the  order  of  0.1  cr2/  r.  The 
computed  diffusion  constant  is  about  one  order  of  magnitude 
higher  than  the  three-bead-chain  model  [11]. 

An  interesting  observation  from  the  phase  diagram  is  that 
the  gel  phase  occurs  at  both  small  and  large  £.  We  attribute 
this  phasic  behavior  to  the  variation  of  the  equilibrium  inter¬ 
particle  distance  d  with  £  due  to  the  second-nearest-neighbor 
effects.  Figure  5(a)  (double  y-axis)  shows  that  d  monotoni- 
cally  increases  with  £,  while  a  maximum  of  D  exists  at  d 
=dc~lcr.  When  d<dc,  the  interparticle  interaction  is  domi¬ 
nated  by  the  repulsive  branch  of  the  potential.  In  this  regime, 
increasing  £  leads  to  a  decrease  in  the  repulsive  force  be¬ 
tween  particles,  giving  rise  to  a  higher  diffusivity.  On  the 
other  hand,  when  d>dc,  the  inter-particle  interaction  is 
dominated  by  the  attractive  branch  of  the  potential.  In  this 
regime,  increasing  £  leads  to  an  increase  in  the  attractive 
force,  giving  rise  to  a  lower  diffusivity.  In  Fig.  5(b),  both  D 
and  d  monotonically  increase  with  temperature  and  exhibit  a 
sudden  rise  at  kBT~0.14e,  implying  a  gel-to-fluid  phase 
transition.  The  thermal  expansion  coefficient  aT  of  the  model 
membrane  is  fitted  at  different  regimes,  as  indicated  in  the 
figure. 

Membrane  tension  X  can  be  calculated  by  applying  the 
virial  formula  to  the  one -particle-thick  fluid  membranes  em¬ 
bedded  in  three-dimensional  (3D)  space  as, 


kJ/S 


FIG.  5.  (Color  online)  Diffusion  constant  D  and  interparticle  distance  d  as  functions  of  the  exponent  £(kBr=0.2e)  (a)  and  temperature  (b) 
for  tensionless  membranes.  The  thermal  expansion  coefficient  aT  of  the  membrane  at  different  temperature  ranges  is  fitted  and  indicated  in 
(b). 
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AA/L2 

FIG.  6.  (Color  online)  Membrane  tension  as  a  function  of  area 
strain  A AIL2.  Area  compression  modulus  KA  is  fitted  to  be  about 
18 k^T/o2.  The  parameters  used  in  the  simulations  are  £=4,  /u=3, 
<90  =  0,  and  kBT=  0.23e. 


2  =  - 


3  NkBT 
2  A 


-2  r,. 
2 A,~  lJ 


(9) 


where  A  is  the  membrane  contour  area,  F(/  is  the  force  ex¬ 
erted  on  particle  i  by  j.  The  first  term  on  the  right-hand  side 
of  Eq.  (9),  denoted  by  2,,  is  the  kinetic  contribution,  while 
the  second  term,  denoted  by  2/,  is  the  potential  contribution. 
The  correctness  of  Eq.  (9)  deserves  some  discussion  here. 
First,  considering  that  a  small  one -particle-thick  membrane 
patch  is  locally  flat,  the  out-of-plane  stress  components  only 
due  to  the  interparticle  forces  are  negligibly  small  because  r,y 
and  Fy  are  almost  in-plane  vectors.  Therefore,  the  coefficient 
of  2 f  is  set  to  1/2.  Second,  for  a  fluid  membrane  in  equilib¬ 
rium,  membrane  tension  is  constant  everywhere.  Therefore, 
the  scalar  tension  calculated  in  Eq.  (9)  around  each  particle 
can  be  summed  and  averaged  over  the  whole  membranes 
even  for  curved  shapes.  Third,  the  coefficient  of  3/2  of  term 
2,  comes  from  the  fact  that  velocity  components  in  all  direc¬ 
tions  contribute  to  expanding  membrane  area.  Figure  6  plots 
the  kinetic  and  potential  contributions  and  the  total  mem¬ 
brane  tension  as  a  function  of  the  imposed  area  strain  ©A 


=  A AIL2.  The  simulations  were  based  on  a  square  membrane 
of  side  length  L~  140cr.  At  a  critical  area  strain  of  —0.09, 
membrane  pore  appears.  This  critical  area  strain  agrees  well 
with  experimental  data  [38].  The  kinetic  and  potential  con¬ 
tributions  to  the  membrane  tension  are  on  the  same  order  of 
magnitude,  suggesting  that  neglecting  the  kinetic  term  would 
lead  to  misinterpretation  of  membrane  tension.  The  resulting 
membrane  tension  is  positive,  and  nearly  linearly  scales  with 
the  imposed  area  strain.  The  slope  of  the  tension-strain  curve 
in  the  large  area  strain  regime  corresponds  to  the  area  com¬ 
pression  modulus,  KA=2/@A.  A  linear  fitting  of  the  curve 
yields  KA~  18kB77 a2. 

In  the  continuum  limit,  the  fluctuation  spectrum  of  a 
square  planar  membrane  of  side  length  L  subjected  to  lateral 
tension  2  takes  the  form  [39], 


kyT 

L\Bq 4  +  2 q2)  ’ 


(10) 


where  B  is  bending  rigidity,  q  is  wave  number.  The  simulated 
fluctuation  spectra  of  a  square  membrane  (consists  of  23595 
particles)  of  side  length  L  —  140cr  at  zero  and  finite  tensions 
were  given  in  Fig.  7(a).  The  q~ 4  dependence  of  the  fluctua¬ 
tion  spectrum  under  zero  tension  is  well  predicted.  At  a  finite 
membrane  tension,  there  exists  a  critical  wave  number  qc 
-sZlB  below  which  the  fluctuation  amplitude  starts  deviat¬ 
ing  from  the  zero-tension  curve  and  exhibits  a  clear  q~ 1  de¬ 
pendence.  Linear  fitting  of  these  two  curves  at  zero  and  finite 
membrane  tensions  yields  the  bending  rigidity  and  mem¬ 
brane  tension,  respectively.  The  fitted  membrane  tension  is 
—  10%  different  from  the  imposed  membrane  tension.  Corre¬ 
sponding  to  the  wavelength  of  — 8o\  there  exists  an  upper 
limit  of  q  =  qm  beyond  which  Eq.  (10)  no  longer  holds.  Both 
qc  and  qm  are  marked  in  Fig.  7(a). 

The  parametric  dependence  of  bending  rigidity  on  //  is 
shown  in  Fig.  7(b),  where  the  bending  rigidity  is  fitted  at 
zero-tension  simulations.  The  bending  rigidity  monotonically 
increases  with  /x.  For  /x  in  the  range  of  2.4-6,  the  obtained 
bending  rigidity  ranges  from  —12 kBT  to  —40 kBT,  which 
falls  in  the  range  of  experimental  data.  There  exists  a  thresh- 


FIG.  7.  (Color  online)  (a)  Fluctuation  spectra  under  zero  and  finite  tensions  for  planar  membranes.  Both  q~ 4  and  q~ 2  dependences  are  well 
predicted  in  simulations.  The  parameters  used  in  the  simulations  are  £=4,  /U=3,  0o  =  O,  and  kBT=  0.23e.  (b)  Membrane  bending  rigidity  (in 
zero-tension  states)  monotonically  increases  with  /x,  while  the  order  of  alignment  of  the  axis  of  symmetry  of  particles  monotonically 
decreases  with  /x.  The  parameters  used  in  the  simulations  are  £=4,  <90  =  0,  and  kBT=  0.23e. 
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FIG.  8.  (Color  online)  Spontaneous  curvature  mediated  phase 
separation  and  subsequent  budding. 

old  value  of  /i  below  which  the  membrane  loses  its  integrity 
due  to  the  weak  orientation  preference.  The  parameter  /i  also 
governs  the  order  of  alignment  of  particle-fixed  vector  n, 
characterized  by  the  average  angle  12  =  cos~l((n, -e,)),  where 
e,  is  the  normal  direction  of  the  reference  plane.  As  shown  in 
Fig.  7(b),  11  ranges  from  about  9°  to  15°,  and  decreases  with 
increasing  bending  rigidity. 

Due  to  the  orientation  dependence  of  the  interparticle  po¬ 
tential,  the  inter-particle  interaction  depends  on  the  direction 
along  which  two  particles  approach  each  other,  which  indi¬ 
cates  that  the  particle  is  geometrically  anisotropic.  The  as¬ 
pect  ratio  can  be  deduced  from  the  bending  and  area  com¬ 
pression  moduli.  For  a  thin  structureless  plate  with  vanishing 
two-dimensional  shear  modulus,  B=KAd2thl  12  holds  [40], 
where  dth  is  the  plate  thickness.  Taking  KA~  1 8AB77 a2  and 
B~20kBT,  the  particle  aspect  ratio  of  this  model  is  found  to 
be  d,hl cr~3.1 .  The  thickness  of  bilayer  membrane  is  dth 
~5  nm,  which  specifies  the  length  scale  of  our  model  cr 

—  1.4  nm  as  far  as  a  bilayer  membrane  is  concerned.  Given 
that  a  single  lipid  molecule  occupies  an  area  of  —0.5  nm2, 
each  particle  in  our  membrane  model  represents  a  few  lipids 
in  the  lateral  direction.  A  typical  value  for  the  diffusion  con¬ 
stant  of  lipids  in  real  phospholipid  membranes  is  about 
1  /im2/s  [41].  The  diffusion  constant  of  our  model  is 
— 0.1  a2 It,  which  maps  out  the  time  scale  of  our  model  t 

—  0.1  /is.  Both  the  length  and  time  scales  are  about  one 
order  of  magnitude  higher  than  the  chain-of-bead  coarse¬ 
grained  model  [11].  Here,  the  length  scale  comparison  is 
based  on  the  fact  that  each  particle  in  our  model  represents  a 
few  lipids  in  two  leaflets  of  the  bilayer.  So  the  degrees  of 
freedom  involved  in  the  same  membrane  area  are  about  one 
order  of  magnitude  fewer  than  the  chain-of-bead  models. 

We  next  demonstrate  the  applicability  of  our  model  by 
selected  examples.  For  bilayer  membranes,  “effective”  spon¬ 
taneous  curvature  may  originate  from  molecule  asymmetry 
[42],  area  mismatch  between  two  leaflets  [43]  of  the  bilayer, 
or  protein-lipid  hydrophobic  mismatch,  or  protein-assisted 
curvatures  [14,44].  Spontaneous  curvature  plays  a  critical 
role  in  determining  shape  transformations  in  vesicles  and  red 


blood  cells  [45,46],  and  is  also  thought  to  be  related  to  lipid¬ 
sorting  in  biological  membranes  [47].  Figure  8  demonstrates 
phase  separation  in  a  binary  lipid  membrane.  The  number 
ratio  of  the  particle  types  A  (red)  and  B  (green)  is  1:10.  We 
assign  different  0Q  and  /i  to  the  different  combinations  of 
particle  pairs,  i.e.,  0$A=  0%B=  1 1.5°,  ^B  =  0°,  /iaa=/iab=6, 
and  /j.bk  =  3.  Other  parameters  used  in  the  simulations  are: 
kBT=  0.23e,  £=4.  It  should  be  noted  that  different  spontane¬ 
ous  curvatures  of  the  two  species  result  in  a  line  tension  that 
induces  the  phase  separation.  Starting  from  random  mixture 
of  the  particle  species,  our  CGMD  simulations  show  that  the 
two  lipid  species  demix  rapidly,  forming  small  domains  (/ 
=  800r).  The  small  domains  gradually  grow  into  stable  do¬ 
mains  of  roughly  the  same  size  (r=4000r).  Further  coales¬ 
cence  of  neighboring  domains  seems  to  be  inhibited  due  to 
the  membrane-mediated  elastic  interaction  [48]. 

Figure  9  shows  four  snapshots  of  our  simulations  of 
adhesion-driven  endocytosis  of  nanoparticles  (NP)  (red).  The 
interaction  between  the  NP  and  the  membrane  particles  is 
nonspecific  and  described  by  the  4-2  LJ  potential  with  a 
cutoff  of  3. Ocr.  The  parameters  for  membrane  particles  are: 
kBT=0.23e,  £=4,  /i=3.0,  and  $0=0°.  The  energy  depth  of 
the  4-2  LJ  potential  is  e,  corresponding  to  —4.4 kBT.  The 
total  resulting  adhesion  energy  (depending  on  the  number  of 
membrane  particles  that  adhere  to  the  NP)  is  sufficient  to 
overcome  the  bending  energy  penalty  (87 tB)  and  drive  NP 
endocytosis  [49-51].  The  membrane  particles  in  blue  are 
within  the  interaction  range  of  the  NP.  For  the  clarity  of 
visualization,  the  front  half  of  the  membrane  model  is  not 
shown.  The  snapshots  show  the  sequence  of  endocytosis  of 
the  NP  (from  left  to  right):  docking,  partial  wrapping,  pinch¬ 
ing  off.  The  NP  is  then  internalized,  while  acquiring  a  layer 
of  membrane  particles. 


IV.  CONCLUSION 

In  conclusion,  we  have  developed  a  particle -based 
solvent-free  fluid  membrane  model.  Several  unique  features 
are  attributed  to  the  success  of  this  model.  First,  the  simple 
mathematical  form  of  the  interparticle  pair  potential  and  the 
one-particle-thick  coarse-graining  substantially  improve  the 
accessible  length  and  time  scales.  The  great  computational 
efficiency  enables  the  model  to  simulate  a  broad  range  of 
membrane  processes  that  occur  at  microscopic  scale.  Second, 
the  pair  potential  is  highly  tunable  and  faithfully  reproduces 
biologically  relevant  membrane  properties.  The  three  model 
parameters  independently  and  effectively  control  diffusion 
constant,  bending  rigidity,  and  spontaneous  curvature,  re¬ 
spectively.  Due  to  the  particle -based  nature,  our  model  simu- 


FIG.  9.  (Color  online)  Simulation  snapshots  of  nanoparticle  endocytosis. 
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lates  molecular  diffusion,  fluidity  and  topological  changes  of 
membranes  naturally. 

Over  past  decades,  considerable  attention  has  been  paid  to 
simulate  vesicle-substrate  interactions  and  interpret  experi¬ 
mentally  measured  membrane  fluctuations  of  human  red 
blood  cells  (RBC)  at  normal  or  disease  states  [52-55].  The¬ 
oretical  analysis  of  membrane  fluctuations  has  been  mainly 
limited  to  planar  membranes  [53]  or  quasispherical  vesicles 
[56].  Triangulated  Monte  Carlo  membrane  models  have  been 
used  to  simulate  fluid  membrane  fluctuations  [23],  where 
“fluidity”  is  realized  by  dynamic  bond  flipping.  The  chain- 
of-bead  models  are  typically  inaccessible  to  long-wavelength 
modes.  As  previously  demonstrated,  our  model  faithfully  re¬ 
produces  fluctuation  spectrum  at  zero  and  finite  tension. 
Therefore,  it  stands  as  an  ideal  model  for  studying  fluctua¬ 
tions  of  membranes  in  gel  or  fluid  phase,  or  membrane  fluc¬ 
tuations  constrained  by  substrate  or  coupled  to  protein  net¬ 
work  [55]. 

Presented  elsewhere  [57],  we  have  established  an  efficient 
algorithm  for  controlling  the  volume  of  a  vesicle  constituted 
by  the  particles.  Reducing  the  vesicle  volume  using  this  al¬ 
gorithm  enables  a  series  of  vesicle  shape  transformations. 


PHYSICAL  REVIEW  E  82,  011905  (2010) 

Combining  our  one-particle-thick  membrane  model  with  the 
volume-control  algorithm  thus  enables  studies  of  fluctuation 
of  non-spherical  vesicles. 

It  has  been  long  unclear  whether  a  one -particle-thick 
membrane  model  with  compatible  membrane  properties  can 
be  established  using  a  simple  anisotropic  pair  potential.  Be¬ 
ing  as  simple  as  LJ  potential,  our  model  for  the  first  time 
provides  a  conclusive  answer.  We  attribute  the  success  of  our 
model  to  two  aspects:  soft-core  repulsive  and  attractive  in¬ 
teractions  and  a  strongly  biased  relative  orientation  between 
a  particle  pair.  The  unique  combination  of  the  simple  math¬ 
ematical  form  and  its  comprehensive  capabilities  in  simulat¬ 
ing  2D  fluid  membrane  physics  can  help  understand  funda¬ 
mental  aspects  of  soft  condensed  matter  physics  of  fluid 
membranes. 

ACKNOWLEDGMENTS 

SLZ  acknowledges  the  support  from  National  Science 
Foundation  under  Awards  No.  0826841  and  No.  0600642.  JL 
acknowledges  the  support  by  NSF  DMR-0520020  and 
AFOSR  FA9550-08-1-0325. 


[1]  L.  F.  Zhang,  K.  Yu,  and  A.  Eisenberg,  Science  272,  1777 
(1996). 

[2]  R.  R.  Gullapalli,  M.  C.  Demirel,  and  P.  J.  Butler,  Phys.  Chem. 
Chem.  Phys.  10,  3548  (2008). 

[3]  M.  Deserno,  Macromol.  Rapid  Commun.  30,  752  (2009). 

[4]  H.  Noguchi,  J.  Phys.  Soc.  Jpn.  78,  041007  (2009). 

[5]  R.  Goetz  and  R.  Lipowsky,  J.  Chem.  Phys.  108,  7397  (1998). 

[6]  M.  Laradji  and  P.  B.  Sunil  Kumar,  Phys.  Rev.  Lett.  93,  198105 
(2004). 

[7]  S.  Yamamoto  and  S.  Hyodo,  J.  Chem.  Phys.  118,  7937  (2003). 

[8]  H.  Noguchi  and  M.  Takasu,  Phys.  Rev.  E  64,  041913  (2001). 

[9]  O.  Farago,  J.  Chem.  Phys.  119,  596  (2003). 

[10]  Z.  I.  Wang  and  D.  Frenkel,  J.  Chem.  Phys.  122,  234711 
(2005). 

[11]  I.  R.  Cooke  and  M.  Deserno,  J.  Chem.  Phys.  123,  224710 
(2005). 

[12]  G.  Brannigan,  P.  F.  Philips,  and  F.  L.  H.  Brown,  Phys.  Rev.  E 
72,  011915  (2005). 

[13]  J.  D.  Revalee,  M.  Laradji,  and  P.  B.  S.  Kumar,  J.  Chem.  Phys. 
128,  035102  (2008). 

[14]  B.  J.  Reynwar  et  al.,  Nature  (London)  447,  461  (2007). 

[15]  I.  R.  Cooke,  K.  Kremer,  and  M.  Deserno,  Phys.  Rev.  E  72, 
011506  (2005). 

[16]  J.  M.  Drouffe,  A.  C.  Maggs,  and  S.  Leibler,  Science  254,  1353 
(1991). 

[17]  G.  Brannigan  and  F.  L.  H.  Brown,  J.  Chem.  Phys.  120,  1059 
(2004). 

[18]  P.  Ballone  and  M.  G.  Del  Popolo,  Phys.  Rev.  E  73,  031404 
(2006). 

[19]  T.  Kohyama,  Physica  A  388,  3334  (2009). 

[20]  H.  Noguchi  and  G.  Gompper,  Proc.  Natl.  Acad.  Sci.  U.S.A. 
102,  14159  (2005). 

[21]  E.  Atilgan  and  S.  X.  Sun,  I.  Chem.  Phys.  126,  095102  (2007). 


[22]  G.  Gompper  and  D.  M.  Kroll,  J.  Phys.:  Condens.  Matter  9, 
8795  (1997). 

[23]  H.  G.  Dobereiner  et  al.,  Phys.  Rev.  Lett.  91,  048301  (2003). 

[24]  J.  Li  et  al.,  Biophys.  J.  88,  3707  (2005). 

[25]  M.  Dao,  J.  Li,  and  S.  Suresh,  Mater.  Sci.  Eng.,  C  26,  1232 
(2006). 

[26]  I.  Li  et  al.,  Proc.  Natl.  Acad.  Sci.  U.S.A.  104,  4937  (2007). 

[27]  F.  Feng  and  W.  S.  Klug,  J.  Comput.  Phys.  220,  394  (2006). 

[28]  H.  Noguchi  and  G.  Gompper,  Phys.  Rev.  E  73,  021903  (2006). 

[29]  L.  T.  Gao,  X.  Q.  Feng,  and  H.  I.  Gao,  J.  Comput.  Phys.  228, 
4162  (2009). 

[30]  T.  Biben,  K.  Kassner,  and  C.  Misbah,  Phys.  Rev.  E  72,  041921 
(2005). 

[31]  Q.  Du,  C.  Liu,  and  X.  Q.  Wang,  J.  Comput.  Phys.  198,  450 
(2004). 

[32]  W.  Humphrey,  A.  Dalke,  and  K.  Schulten,  J.  Mol.  Graphics 
14,  33  (1996). 

[33]  J.  G.  Gay  and  B.  J.  Berne,  J.  Chem.  Phys.  74,  3316  (1981). 

[34]  R.  Everaers  and  M.  R.  Ejtehadi,  Phys.  Rev.  E  67,  041710 
(2003). 

[35]  S.  Nose,  Mol.  Phys.  52,  255  (1984). 

[36]  W.  G.  Hoover,  Phys.  Rev.  A  31,  1695  (1985). 

[37]  H.  J.  C.  Berendsen  et  al.,  J.  Chem.  Phys.  81,  3684  (1984). 

[38]  R.  Lipowsky  and  E.  Sackmann,  Structure  and  Dynamics  of 
Membranes  (Elsevier  Science,  Amesterdam,  1995),  Vol.  1  and 
2. 

[39]  F.  Brochard  and  J.  F.  Lennon,  J.  Phys.  (Paris)  36,  1035  (1975). 

[40]  R.  Goetz,  G.  Gompper,  and  R.  Lipowsky,  Phys.  Rev.  Lett.  82, 
221  (1999). 

[41]  P.  F.  Fahey  and  W.  W.  Webb,  Biochemistry  17,  3046  (1978). 

[42]  W.  Helfrich,  Z.  Naturforsch  C  28,  693  (1973). 

[43]  M.  P.  Sheetz  and  S.  J.  Singer,  Proc.  Natl.  Acad.  Sci.  U.S.A.  71, 
4457  (1974). 


011905-7 


YUAN  et  al. 


PHYSICAL  REVIEW  E  82,  011905  (2010) 


[44]  P.  D.  Blood  and  G.  A.  Voth,  Proc.  Natl.  Acad.  Sci.  U.S.A.  103, 
15068  (2006). 

[45]  H.  W.  G.  Lim,  M.  Wortis,  and  R.  Mukhopadhyay,  Proc.  Natl. 
Acad.  Sci.  U.S.A.  99,  16766  (2002). 

[46]  U.  Seifert,  Adv.  Phys.  46,  13  (1997). 

[47]  U.  Seifert,  Phys.  Rev.  Lett.  70,  1335  (1993). 

[48]  T.  S.  Ursell,  W.  S.  Klug,  and  R.  Phillips,  Proc.  Natl.  Acad.  Sci. 
U.S.A.  106,  13301  (2009). 

[49]  H.  Y.  Yuan  and  S.  L.  Zhang,  Appl.  Phys.  Lett.  96,  033704 
(2010). 

[50]  S.  L.  Zhang  et  al.,  Adv.  Mater.  21,  419  (2009). 


[51]  M.  Deserno  and  T.  Bickel,  EPL  62,  767  (2003). 

[52]  H.  Strey  and  M.  Peterson,  Biophys.  J.  69,  478  (1995). 

[53]  T.  Auth,  S.  A.  Safran,  and  N.  S.  Gov,  Phys.  Rev.  E  76,  051910 
(2007). 

[54]  Y.  Park  et  al. ,  Proc.  Natl.  Acad.  Sci.  U.S.A.  107,  1289  (2010). 

[55]  Y.  K.  Park  et  al.,  Proc.  Natl.  Acad.  Sci.  U.S.A.  105,  13730 
(2008). 

[56]  S.  T.  Milner  and  S.  A.  Safran,  Phys.  Rev.  A  36,  4371  (1987). 

[57]  H.  Yuan,  C.  Huang,  and  S.  Zhang,  Soft  Matter  (to  be  pub¬ 
lished). 


011905-8 


